clear all;
load synthpre.mat
load lookupt.mat

% Compute the table of optimal reserve price for each std/mean ratio that we consider

RPopt = [];
m = 1;

for stdrat = 0.1:0.1:1
    v = (m*stdrat)^2;
    mu = log((m^2)/sqrt(v+m^2));
    sigma = sqrt(log(v/(m^2)+1));
    fun = @(x)(-max(0.00001,x)*(1-normcdf(log(max(0.00001,x)),mu,sigma)));
    x0 = 1;
    rstar = fminsearch(fun,x0);
    RPopt = [RPopt; rstar];
end

estrp = zeros(size(synthpre,1),1);

for kw = 1:size(synthpre,1)

    % The three observed moments here

    nbids = synthpre(kw,8);
    meanbid = synthpre(kw,9);
    stdbid = synthpre(kw,10);
    
    % Match the closest profile of primitives
    
    j = find(lookupt(:,1) >= nbids, 1);
    n = size(lookupt,1);
    
    % Weighted sum of squared deviations

    delta = 0.01*(nbids - lookupt(j:n,4)).^2 + (meanbid - lookupt(j:n,5)).^2 + (stdbid - lookupt(j:n,6)).^2;
    
    [v,i] = min(delta);
    
    lookupt(i+j-1,:);
    
    eststdrat = round(10*lookupt(i+j-1,3));
    
    rstar = lookupt(i+j-1,2)*RPopt(eststdrat);

    estrp(kw,:) = [rstar];

end

estrp = [synthpre estrp];

% Save the estrp table, both as a text file and a Matlab file
writematrix(estrp);
save estrp.mat estrp

% Output the data corresponding to Figure 1 and Table 4
edges = 0:0.05:2.5;
histogram(estrp(:,size(estrp,2)),edges,'Normalization','probability')

[prctile(estrp(:,size(estrp,2)),10) prctile(estrp(:,size(estrp,2)),25) prctile(estrp(:,size(estrp,2)),50) prctile(estrp(:,size(estrp,2)),75) prctile(estrp(:,size(estrp,2)),90)]
